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We use molecular dynamics simulations for a first principles-based effective Hamiltonian to cal¬ 
culate two important quantities characterizing the electrocaloric effect in BaTiC> 3 , the adiabatic 
temperature change AT and the isothermal entropy change AS, for different electric field strengths. 

We compare direct and indirect methods to obtain AT and AS, and we confirm that both methods 
indeed lead to identical result provided that the system does not actually undergo a first order phase 
transition. We also show that a large electrocaloric response is obtained for electric fields beyond 
the critical field strength for the first order phase transition. Furthermore, our work fills several 
gaps regarding the application of the first principles-based effective Hamiltonian approach, which 
represents a very attractive and powerful method for the quantitative prediction of electrocaloric 
properties. In particular, we discuss the importance of maintaining thermal equilibrium during the 
field ramping when calculating AT using the direct method within a molecular dynamics approach. 


I. INTRODUCTION 


The ongoing search for alternative cooling technolo¬ 
gies which are more energy-efficient and environmentally 
friendly than conventional vapor-compression refrigera¬ 
tors and offer the additional possibility for device minia¬ 
turization has boosted research activities within the fields 
of electrocaloric, elastocaloric, and magnetocaloric ef- 
fectsii~— The common feature in all three cases is that 
the application of an external field (either electric, stress, 
or magnetic field) under adiabatic conditions, i.e. when 
the active material is thermally isolated from the envi¬ 
ronment, results in a temperature change of the corre¬ 
sponding material. This reversible temperature change 
can be used to transfer heat from a cool reservoir (the 
heat load) to a warmer reservoir (e.g. the environment), 
thereby lowering the temperature of the heat load (or 
keeping it at constant low temperature). It has been 
found that such caloric effects are especially large close to 
ferroic first order phase transitions, where giant responses 
can be triggered through relatively modest fields^ - — 

In particular the electrocaloric (EC) effect has become 
very attractive for potential future applications, due to 
the discovery of a giant EC temperature change of 12 K 
in Pb(Zr,Ti )03 thin fflms4 Here, the high crystalline 
quality that can be achieved in thin film samples allows 
the application of rather high electric fields without trig¬ 
gering a dielectric breakdown of the samples. In recent 
years, a large number of studies - both theoretical and 
experimental have contributed to a better understand¬ 
ing of the EC effect (see, e.g. Refs. and references 
therein). 

Nevertheless, direct measurements of the adiabatic 
temperature change are still rather challenging, in par¬ 
ticular for the case of thin film samples. Therefore, an 
indirect determination of this temperature change is of¬ 
ten preferred. The indirect method is based on a ther¬ 


modynamic Maxwell relation connecting the isothermal 
field-induced entropy change with the temperature de¬ 
pendence of the electric polarization at fixed electric field: 



T 



( 1 ) 


The adiabatic temperature change AT can then be ob¬ 
tained from pyroelectric measurements, i.e. by measuring 
the electric polarization P as function of temperature T 
at different electric fields £: 


AT = - 





( 2 ) 


Here, C p £ is the specific heat at constant pressure and 
applied field, and the external field is varied from £\ to £ 2 - 
It has to be noted that, if the system undergoes a first or¬ 
der phase transition, the derivative dP/dT is ill-defined 
and the specific heat diverges, which in principle does 
not allow application of Eq. ©• Furthermore, a possible 
contribution to the EC effect stemming from the latent 
heat of the first order phase transition is not accounted 
for by Eq. ©• Instead, the Clausius-Clapeyron equation 
has to be used to obtain the corresponding contribution. 
In addition, the indirect method is only suitable for er- 
godic systems. For example, it was shown that the results 
from direct and indirect measurements do not match for 
relaxor polymers^ but compare well for “normal” ferro¬ 
electric polymers *2 Finally, the influence of domains and 
anisotropy effects are not covered by the scalar form of 
the Maxwell relation, Eq. 

Another important quantity for characterizing the EC 
effect is the isothermal entropy change AS, which is re¬ 
lated to the amount of heat that is required to keep the 
system at constant temperature while an electric field is 
applied or removed. The isothermal entropy change can 
also be obtained indirectly from pyroelectric measure- 
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ments by simply integrating Eq. 0: 


AS = 




( 3 ) 


On the other hand, AS can also be obtained in a (quasi-) 
direct way from integrating the specific heat at constant 
electric field: 

AS = [ Cp ’ £l ~ / Cp ’ S2 d T' . (4) 

J Ti T 

Note that strictly speaking this relation is only valid for 
Ti —> 0. Nevertheless, for sufficiently low Ti one can 
assume that S(T\,£\) ~ S(Ti,£ 2 ) and then Eq. © can 
be expected to give a good estimate of AS 

In the work presented in this article, we use a first 
principles-based effective Hamiltonian approach^— to 
calculate the EC effect in the prototypical ferroelectric 
perovskite BaTiC> 3 , and to address the applicability of 
the indirect method for evaluating AT and AS. Perform¬ 
ing micro-canonical molecular dynamics (MD) on the ef¬ 
fective Hamiltonian allows for a direct calculation of the 
EC temperature change under application or removal of 
an electric field. Within the same framework, the temper¬ 
ature dependence of the electric polarization under differ¬ 
ent electric fields can be calculated and the temperature 
and entropy changes can then be evaluated via Eqs. © 
and ©■ Thus, the effective Hamiltonian provides a sim¬ 
plified but nevertheless realistic “testing ground” for the 
general applicability of the indirect methods. 

Previous studies employing first principles-based effec¬ 
tive Hamiltonians have found good agreement between 
direct and indirect calculations of the EC temperature 
change ) 15 ' 16 provided that both £\ and £ 2 are above the 
critical field for the first order phase transition, i.e. in a 
regime where no discontinuities of the polarization oc¬ 
cur as function of temperature and electric field. In 
Ref. [H the EC temperature change for Ba 0 . 5 Sr 0 . 5 TiO 3 
has been calculated using micro-canonical Monte Carlo 
simulations (Creutz algorithm), and the so-obtained val¬ 
ues have been compared with the indirect evaluation 
based on Eq. ©, where P(T,£) has been obtained from 
standard Monte Carlo simulations within the canonical 
ensemble. In Ref. El, the direct calculation of AT for 
BaTiC >3 has been performed using a micro-canonical MD 
algorithm. The indirect evaluation of AT using Eq. © 
showed reasonable agreement with the corresponding di¬ 
rectly calculated values. Discrepancies were attributed to 
inconsistencies arising from an empirical, temperature- 
dependent, pressure correction and to the use of a con¬ 
stant empirical value for the specific heat (the experi¬ 
mental value for C p £ at room temperature was used in 
Ref. (TI). It is important to note that C Pt £ is not con¬ 
stant and varies significantly with temperature and ap¬ 
plied field, especially near the phase transition* 17 

Here, we calculate the specific heat of the effective 
Hamiltonian, as function of temperature and electric 
field, in order to allow for a fully consistent compari¬ 
son between the direct and indirect evaluation of AT 


and AS. We confirm that both methods indeed lead to 
identical result provided that the system does not ac¬ 
tually undergo a first order phase transition. We also 
show that the actual transition is not crucial for obtain¬ 
ing a sizable EC response and compare this with the case 
of magnetocaloric Heusler alloys. Furthermore, we calcu¬ 
late the isothermal EC entropy change and again demon¬ 
strate good agreement between direct and indirect meth¬ 
ods. Finally, we investigate how fast the electric field can 
be changed within the MD simulation without the sys¬ 
tem going out of thermal equilibrium. In particular, we 
demonstrate the importance of maintaining equilibrium 
during the simulation by monitoring changes in the total 
energy of the system. 

This paper is organized as follows. In Section I© we 
briefly describe our computational method. Our results 
are then presented in Section mu which is divided into 
two parts, the first describing the effect of different ramp¬ 
ing rates for the electric field, the second discussing the 
EC temperature and entropy changes. Finally, in the last 
section, we summarize our main results and conclusions. 


II. COMPUTATIONAL METHOD 

For our study, we use the effective Hamiltonian pro¬ 
posed by Zhong ef al. j 12 i 13 This effective Hamiltonian is 
applicable to ferroelectrics with a cubic perovskite parent 
structure. The ferroelectric polarization in these materi¬ 
als can be described by a relative displacement of cations 
and anions, represented by a soft mode variable in the 
Hamiltonian. In addition, local strain variables are in¬ 
cluded. This type of description retains the dominant 
terms in the total energy while reducing the number of 
degrees of freedom per unit cell from 15 to 6 (3 soft mode 
variables and 3 local strain variables). 

All parameters for the effective Hamiltonian can be 
obtained using ab initio density functional theory calcu¬ 
lations*^^ The effective Hamiltonian approach is there¬ 
fore able to determine temperature-dependent properties 
of ferroelectric materials without the need for empirical 
input parameters. For example, it was demonstrated that 
the three consecutive phase transitions in bulk BaTiOs 
are successfully reproduced^ Furthermore, the effective 
Hamiltonian approach has been used successfully for the 
calculation of EC properties, 15 i 16 i 19 ~ — 

We perform MD simulations employing the ef¬ 
fective Hamiltonian as implemented in the feram 
codoi^ (http://loto.sourceforge.net/feram/), using 
the available parameter set for BaTiOsr^ which has been 
obtained using the generalized gradient approximation 
for the exchange-correlation functional according to Wu 
and Cohen^ 

In order to directly calculate the adiabatic EC temper¬ 
ature change, we first thermalize the system at a given 
temperature and electric field using a Nose-Poincare 
thermostat^ We then switch off the thermostat, i.e. 
we switch to the micro-canonical ensemble, and slowly 
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change the electric field while monitoring the resulting 
changes in the total and kinetic energies. These calcula¬ 
tions are performed using a 96 x 96 x 96 supercell, i.e. 
corresponding to 96 simple perovskite unit cells along 
each cartesian direction. A time step of 1 fs per MD step 
is used and the thermalization (averaging) time for these 
direct calculations is equal to 80 ps (40 ps). As usual, the 
temperature is calculated from the kinetic energy, E k j n , 
of the system: 


2A kin 

Nfks 


( 5 ) 


where Nf denotes the number of degrees of freedom of 
the system and ks is the Boltzmann constant. The EC 
temperature change AT is then simply obtained from the 
difference between the initial and final temperature of the 
system, i.e. before and after the electric field is ramped 
on or off. 

To reduce the computational effort, we use a simpli¬ 
fied treatment for the local strain variables, which are 
obtained by minimization of the total energy for the cur¬ 
rent soft-mode configuration in each MD step. Thus, our 
model contains only Nf = 3 dynamic degrees of free¬ 
dom per unit cell (the 3 soft mode variables), compared 
to the original 15. As a result, the model specific heat 
and the directly calculated AT need to be rescaled before 
comparing to experimental dataJ£ However, since the fo¬ 
cus of this work is on the internal consistency within the 
model description, in order to assess the general validity 
of the indirect determination of AT and AS, and not on 
a quantitative comparison with experimental data, we do 
not perform such rescaling within this work, i.e. except 
where otherwise noted, all presented values for AT and 
Cp,£ refer to the model system and not to the real mate¬ 
rial. We also note that a simple rescaling of AT neglects 
the fact that the simulation corresponds to a “wrong” 
final state of the system, i.e. with different temperature 
and polarization compared to the final state that would 
be obtained in a real experiment, and therefore corrects 
only partially for the missing degrees of freedom. 

To calculate the adiabatic temperature change and the 
isothermal entropy change using the indirect method, 
we calculate polarization as function of temperature (on 
a 1K grid) at several applied electric fields using a 
16 x 16 x 16 simulation cell, a thermalization time of 
120 ps and an averaging time of 80 ps, with a 2fs time 
step per MD iteration. These calculations are performed 
in the canonical ensemble using the Nose-Poincare ther¬ 
mostat. We then use smoothing cubic spline functions to 
fit the polarization versus temperature data, in order to 
determine ( 8P/dT) £ . 

The specific heat of the model Hamiltonian at constant 
pressure and electric field, required for the indirect cal¬ 
culation of AT and the (quasi-) direct calculation of AS, 
is determined by calculating the derivative of the total 
energy, i.e. by using the relation C Pi £ = (dE tot /dT) p£ , 
which is applicable for our simulations performed at zero 
pressure. To calculate E tot (T), we use a 96x96x96 sim¬ 


ulation cell, equilibration and averaging times of 80 ps 
and 40 ps, respectively, and a 2fs time step. The tem¬ 
perature dependence of C p> £ has been calculated using 
“cooling” as well as “heating” simulations, i.e. where the 
system at a particular temperature is initialized from a 
thermalized configuration at sightly higher or lower tem¬ 
perature, respectively (see, e.g. Ref. [HI). While an ap¬ 
preciable thermal hysteresis is obtained for zero electric 
field, the thermal hysteresis completely vanishes for fields 
above 20-30 kV/cm. Therefore, only results from “cool¬ 
ing” runs are presented in the following. In the vicin¬ 
ity of the phase transition, due to the sharp features in 
Cp £, a dense IK mesh and extended equilibration time 
is used for field strengths below 75 kV/cm. Otherwise, 
a temperature grid of 5 K is used and the specific heat 
is extrapolated to a 1K temperature grid and a moving 
average is used to further smooth the data. Above 450 K 
and for field strengths of more than 200 kV/cm, the to¬ 
tal energy varies only weakly. Therefore, we have used a 
coarser temperature grid of 10 K in that region. 

Using ( dP/dT) e and the calculated C Pi £, we can then 
obtain AT from Eq. (J2]) and AS from Eq. (J4J) . We 
note that we have confirmed the absence of noticeable 
finite size effect in our results for electric fields above 
~25 kV/cm, which is above the critical field for the first 
order phase transition. Therefore, using ( dP/8T) £ and 
U Pi £ obtained from different sizes of the simulation cell 
does not introduce any significant errors or inconsisten¬ 
cies to our analysis. 

In addition, we have performed test calculations as¬ 
sessing the effect of different field ramping rates (see 
Sec. HITAl) . These tests are performed using a 48 x 48 x 48 
simulation cell and a time step of 1 fs. Different thermal¬ 
ization and averaging times have been used in these cal¬ 
culations, depending on the specific field strength, ramp¬ 
ing rate, and temperature. In all cases we verified that 
the system is sufficiently equilibrated and averages were 
obtained with good accuracy. 

We note that in our calculations, we do not apply any 
empirical pressure corrections, which have been used in 
previous studies to correct for deficiencies of the first 
principles calculations or to mimic thermal expansion. 
Such pressure corrections can lead to better agreement 
between the calculated and measured transition temper¬ 
atures^ However, as already pointed out in Ref. [l6j, a 
temperature-dependent pressure correction can also lead 
to inconsistencies between the direct and indirect calcu¬ 
lation of AT. Consequently, we refrain from using such 
pressure corrections (or from rescaling the parameter H 2 
in the soft mode energy of the effective Hamiltonian, see 
e.g. Ref M) in this work. As a result, our calculated tran¬ 
sition temperature T c for the cubic to tetragonal phase 
transition (~ 270 K) deviates from the known experimen¬ 
tal value (403 K). However, it can be expected that nev¬ 
ertheless trends are accurately described and that the 
calculated temperature changes (after rescaling for the 
correct number of degrees of freedom) are also quantita¬ 
tively of the right magnitude. 
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Case 


Switching on 


Switching off 


Paraelectric phase (T > T c ) 


Instantaneous 0 
Ramping ~^x£ 


2 

app 


X^"app 

±y £ 2 

o X'-'a 


Ferroelectric phase (T < T c ) 
Instantaneous —Po£ apP Po£ apP 

Ramping -P 0 £ app - \x'£l PP Po£ apP + \x! £ 


■ + x!£ 


2 

app 


2 

app 


TABLE I. The changes in the total energy on varying the ap¬ 
plied electric field are tabulated for instantaneous switching 
and slow ramping of the field. “Switching on” corresponds 
to the field varying from zero to £ app , and vice-versa for the 
“switching off” case. The formulas are derived using the fol¬ 
lowing simplified assumption: the induced polarization Pi n d 
depends linearly on the applied field £, the proportionality 
constant is the dielectric susceptibility \ ' n the paraelectric 
phase and x i n the ferroelectric phase. Po is the spontaneous 
polarization of the system in the ferroelectric phase. 


III. RESULTS AND DISCUSSION 
A. Rate dependence 

First, we investigate the influence of the rate of change, 
d£ / dt , with which the electric field is ramped up or down 
in our simulations. This is an important technical point, 
since, depending of course on the invested computational 
resources, MD simulations can only cover time periods 
of up to a few nano seconds. This means that within 
the simulations, the electric held needs to be changed 
extremely fast compared to a real experiment. Neverthe¬ 
less, it is very important to ensure that the system always 
stays in thermal equilibrium and that the MD simulation 
indeed describes a reversible process. 

We start by analyzing the change of the total energy 
under application of an electric held for the two cases 
of instantaneous electric held switching and very slow 
ramping. In general, the change in total energy AP t ot 
under application or removal of an electric held is given 
by AE tot = — f P ■ d£ 1 where P is the polarization of 
the system. For instantaneous switching, the polariza¬ 
tion cannot follow the change of the applied held and 
stays essentially constant during the switching process. 
In the paraelectric phase, the spontaneous polarization 
is zero. Therefore, when the held is switched on in¬ 
stantaneously for T > T c , AE to t is also equal to zero. 
However, if the held is instantaneously switched off from 
some hnite value £ app , then even at T > T c , there is an 
induced polarization, Pj n d = x^app, resulting in a non¬ 
zero AE tot = Find • £ app = X'^app- This implies that the 
complete cycle of applying and removing an electric held 
instantaneously to the system results in an irreversible 
process. 

On the other hand, if the held is applied/removed 
slowly, then the polarization can follow the external held, 
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FIG. 1. (Color online) (a) Schematic depiction of the sim¬ 
ulation cycle (see text). Panels (b) and (c) show the total 
energy as a function of MD steps for instantaneous field appli¬ 
cation/removal and for slow held ramping, respectively. Here, 
the starting temperature Ti is 530 K and the applied held is 
100 kV/cm. In (c) the rate of change of the applied held is 
equal to 0.05 kVcrrUHs -1 . For clarity, the total energy is plot¬ 
ted only after the system is switched to the microcanonical 
ensemble. 


and at each time P = \£. The resulting change in total 
energy is then given according to AE to t = — J £ . P-d£ = 
±|x^a PP - Here, £\ and £f are the initial and hnal applied 
helds, respectively, which are equal to zero and £ app for 
application of the held, and the other way round for re¬ 
moval. The plus and minus signs then correspond to 
removal and application of £ a pp, respectively. Thus, it 
can be seen that slow ramping results in the same mag¬ 
nitude of the total energy change for switching the held 
on and off, i.e. one obtains a reversible process. Similar 
arguments hold true within the ferroelectric phase, with 
an additional term coming from the spontaneous polar¬ 
ization Po. The resulting total energy changes for the 
various cases are tabulated in Table Q] 

Next, we perform simulations at different temperatures 
to examine whether the simple considerations outlined in 
the preceding paragraphs are consistent with the actual 
MD simulations for the effective Hamiltonian. We have 
selected two temperatures, T = 530 K (in the paraelec¬ 
tric phase) and T = 270 K (in the ferroelectric phase, just 
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FIG. 2. (Color online) The change in the total energy 
on switching off the applied field is plotted as a function 
of applied electric field for (a) the paraelectric phase at 
Ti = 530 K and (b) the ferroelectric phase at T) = 270 K. 
The rate of change of the applied field d£/dt is equal to 
—0.05 kVcm _1 fs _1 for the “ramping” case (red squares). The 
symbols show the data points, whereas the dashed lines are 
fits to the corresponding functional forms listed in Table Q] 


below the transition temperature), at which we monitor 
the change in the total energy on switching the field on 
and then off again for several values of £ app . The full 
simulation cycle is depicted schematically in Fig. Oja). 
First, the system is thermalized at temperature 7] and 
field £ = 0 within the canonical ensemble. The simu¬ 
lation is then switched to the microcanonical ensemble, 
the electric field is ramped up to £ app , and the resulting 
temperature change A T on is monitored. Then, the field is 
ramped down again, and the corresponding temperature 
change AT 0 g is monitored. If the system stays in thermal 
equilibrium during the entire simulation cycle, then both 
its total energy and its temperature, Tf, at the end of 
the simulation should be identical to the corresponding 
starting values, and A T a g = — AT on . 

The evolution of the total energy over a full cycle at 
T{ = 530 K, i.e. in the paraelectric phase, is shown in 
Fig. mb). In this simulation, the field is switched in¬ 
stantaneously. As expected, there is no change in the 
total energy while switching on the field (see Table |U), 
but there is a jump of the total energy when the held is 
switched off. Fig. [l)c) shows the evolution of the total 
energy when the held is ramped up and down slowly. In 
this case, |Ai?tot| is the same for application and removal 
of the held. This confirms that very fast switching of the 
applied held results in an irreversible process. 

Further, we plot the change in the total energy AE tot 


as a function of the applied held £ app at Ti = 530 K 
and T[ = 270 K in Figs. EJa) and (b), respectively. 
These AE to t values correspond to removal of the held 
(“switching off”). For slow ramping, these are equal 
to those obtained from switching on (but with oppo¬ 
site sign). The data from the simulations is htted using 
the corresponding functional forms given in Table |U The 
hts are indicated by dashed lines and match very well 
with the data. In the paraelectric phase the total en¬ 
ergy change for switching off the held depends quadrat- 
ically on the applied held strength and there is a fac¬ 
tor of 2 difference between slow ramping and instanta¬ 
neous switching. The ht to the instantaneous switch¬ 
ing data gives x = 3.5 x 10 -2 fj ,C ■ kV^ 1 cm _1 which 
matches well with the corresponding value obtained from 
the ramping data (y = 3.8 x 10~ 2 fiC ■ kV^cm -1 ). In 
the ferroelectric phase, AE tot is dominated by the lin¬ 
ear contribution stemming from the spontaneous polar¬ 
ization. From the ht of the instantaneous switching 
data, we obtain P 0 = 30.32 fiC/cm 2 and x' = 1-3 x 
10~ 2 /xC-kV _1 cm _1 , whereas for the case of slow ramp¬ 
ing the corresponding quantities are 30.32 ^tC/cm 2 and 
1.5 x 10 -2 /iC-kV _1 cm _1 respectively. There is good 
agreement between these two data sets. Similarly, the 
value for Pq obtained from htting AE tot for instanta¬ 
neous “switching on” of the electric held (not shown here) 
is equal to 30.31/xC/cm 2 . We can also compare these 
values for the spontaneous polarization to that obtained 
directly from the MD simulations at T = 270 K, which 
is equal to 28.2 ^C/cm 2 . Note that the agreement be¬ 
tween the various parameters is excellent considering the 
simplicity of the approach. This shows that the simple 
considerations outlined at the beginning of this section 
do indeed lead to a consistent description of the various 
switching cases. 

Fig. (31(a) shows the obtained temperature changes 
when the electric field is switched on and off, AT on and 
AT 0 ff, as function of the inverse rate (d£ / dt) , with 
which the field is ramped up and down, for a start¬ 
ing temperature T) = 350 K (i.e. close to the maxi¬ 
mum EC effect). Instantaneous switching corresponds 
to (dS/dt)- 1 = 0 and slow ramping corresponds to a 
large inverse rate. It can be seen that for small inverse 
ramping rates, i.e. for fast switching, |AT on | ^ |AT 0 ff|- 
The difference between the values is about 9K which 
is significant compared to the converged EC tempera¬ 
ture change of 17.6 K observed at this temperature^ 
This implies that the system goes out of equilibrium 
during fast switching, consistent with the total energy 
considerations discussed above. As the rate is reduced 
(i.e. the inverse rate is increased), the difference be¬ 
tween |AT on | and |AT 0 ff| becomes smaller, and already 
for an inverse rate of 500 kV _1 -fs-cm (corresponding 
to \d£/dt\ = 0.002kVcm _1 fs _1 ), the difference becomes 
negligible. Similar behavior can be observed also for 
other initial temperatures. 

In Fig. (3Kb), the directly calculated EC temperature 
change obtained by instantaneously switching off the 
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FIG. 3. (Color online) (a) EC temperature changes for switch¬ 
ing the electric field on and off, AT on and AT/g, as func¬ 
tion of the inverse rate of change of the applied electric field, 
( d£/dt)~ , for a starting temperature T\ = 350K. The ap¬ 
plied field £ app is equal to 200 kV/cm. The simulation cycle 
used to obtain AT on / off is depicted in Fig. HJa). (b) The EC 
temperature change AT as a function of temperature is plot¬ 
ted for slow ramping and instantaneous ramping of the field. 
The field is varied from 225kV/cm to OkV/cm. 


electric field is compared with the one obtained using 
slow field ramping (d£/dt = —0.002 kV • cm _1 fs _1 ) for 
an initially applied field of 225kV/cm. Note that for 
these calculations the system is thermalized in the canon¬ 
ical ensemble with a nonzero field, which is removed af¬ 
ter switching to the microcanonical ensemble. While the 
overall behavior is the same for both cases, the instanta¬ 
neous switching underestimates the magnitude of AT at 
essentially all temperatures. This is consistent with our 
earlier observations in Fig. ©a). In addition, the temper¬ 
ature for which the largest temperature change occurs is 
slightly shifted to lower temperatures. 

The maximum AT, observed at 320 K in Fig. (3Kb), 
is about — 25.6 K for slow ramping. After appro¬ 
priate rescaling for the correct number of degrees of 
freedom. ) 16 ! 26 this corresponds to an EC temperature 
change of around 5.1 K, which agrees well with earlier 
reports for similar applied field strengths.™^ After scal¬ 
ing, the difference between AT calculated using instanta¬ 
neous switching and using slow field ramping is approxi¬ 
mately 1K for the given field strength (except very close 
to the peaks). 

Our results up to now thus demonstrate the necessity 
of ensuring that the system is in thermal equilibrium 
throughout the whole MD simulation for a correct direct 
calculation of the adiabatic EC temperature change. In 
the following we use a rate d£/dt = 0.002 kVcm _1 fs _1 for 
all our direct calculations of the EC temperature change. 
Even though this rate is very fast compared to actual 
experimental rates, it is sufficiently slow to avoid irre¬ 
versibility in the calculation and also allows for reason¬ 
able simulation times. 



T (K) 


FIG. 4. (Color online) Calculated specific heat of the model 
Hamiltonian as a function of temperature for different applied 
electric fields. For better comparison, the Dulong-Petit value 
of 3 fcs is indicated by the thick horizontal black line and 
the obtained high temperature limit is indicated by the thin 
horizontal red line. 

B. Direct versus indirect EC effect 

Next, we compare results from the direct and indirect 
approaches. As pointed out previously, all calculations 
are performed using the same effective model Hamilto¬ 
nian and no experimental data is used. In our previ¬ 
ous work, we have calculated the EC effect using the 
indirect method,— 1 ' 22 but we have used the experimen¬ 
tal specific heat value at room temperature to evaluate 
Eq. ©. Although this is a valid first approximation, this 
treatment ignores the temperature and electric field de¬ 
pendence of C p ^s , as well as the mismatch between the 
number of degrees of freedom of the real system and the 
model, which will lead to differences between the results 
obtained from direct and indirect methods. Therefore, 
in the present work, we calculate the specific heat from 
the model Hamiltonian as a function of temperature at 
different applied fields. This allows for an internally con¬ 
sistent comparison between direct and indirect methods, 
and also enables us to obtain A S dlr using the (quasi-) 
direct method, Eq. 

In Fig. [4] we show our results for the specific heat 
of the effective Hamiltonian as a function of tempera¬ 
ture at several applied fields. In absence of an applied 
field, there is a pronounced peak (divergence) at the fer¬ 
roelectric transition, which also shows pronounced ther¬ 
mal hysteresis. Such a divergence is characteristic for 
a first order phase transition. With increasing electric 
field, the phase transition and thus the peak in the spe¬ 
cific heat shift to higher temperature. Furthermore, the 
transition becomes smoother and the thermal hysteresis 
disappears for fields of around 25kV/cm and stronger. 
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FIG. 5. (Color online) (a) Comparison of the EC tempera¬ 
ture change as function of the initial temperature obtained 
using direct and indirect methods. The field is varied from 
300kV/cm to 75kV/cm. In the direct calculation slow field 
ramping with dS/dt = —0.002 kVcm _1 fs _1 is used, (b) The 
EC entropy change as a function of temperature obtained 
from direct and indirect methods. Here, the applied field is 
varied from 75kV/cm to 300kV/cm. The results of the di¬ 
rect calculations, but with an offset such that it matches the 
indirect calculation at T = 200 K (see text) are indicated by 
the blue dotted line. 


Previous phenomenological thermodynamic calculations 
and experimental measurements on bulk BaTiC >3 have 
shown that an applied field of about lOkV/cm is suffi¬ 
cient to suppress the first order phase transition. 1 '' 2 ' 28 

Both at high and low temperatures, i.e. away from the 
phase transition, C Pi s approaches constant values. We 
recall that our simulations contain only 3 degrees of free¬ 
dom instead of 15 for the real system. If the Hamiltonian 
would be exactly quadratic in the 3 soft mode variables, 
then each degree of freedom would contribute Iks to the 
specific heat, and we would expect a value of 3 ks, which 
is equal to 0.65 J-K _1 cnC 3 , both in the high and low 
temperature limit. We note that in our purely classical 
simulations no modes are “frozen in” at low tempera¬ 
tures. We find a high-temperature limit of 0.59 J/K/cm 3 
from our calculations. This small discrepancy with the 
Dulong-Petit value of 3 kg results from the higher order 
terms in the Hamiltonian. These contribute less at low 
temperatures, and therefore at low temperatures below 
T c the calculated value compares well with the Dulong- 
Petit value. 

Next, we compare both the calculated adiabatic EC 
temperature change AT as well as the isothermal EC 
entropy change AS obtained using direct and indi¬ 
rect methods. As mentioned before, we use a rate of 
—0.002 kVcm _1 fs _1 for changing the applied field. In 
Fig. [U)a), we compare the results from direct and indi¬ 
rect calculations of the adiabatic EC temperature change 
AT. These results correspond to removal of the field 
(“switching off”), i.e. a negative AT (as in Fig. [3|b)). 



FIG. 6. (Color online) EC temperature change AT as func¬ 
tion of temperature for different field intervals with the same 
total width of 50kV/cm. Results obtained using the direct 
(indirect) method are shown as dotted lines with symbols 
(solid lines without symbols). 


We consider only fields > 75kV/cm for this comparison, 
in order to exclude the region close to the first order 
phase transition, where the indirect method is not appli¬ 
cable. It can be seen that the AT values calculated from 
direct and indirect methods match extremely well over 
the whole temperature range. This shows that within 
our consistent description, where all quantities are cal¬ 
culated using the same effective Hamiltonian, the direct 
and indirect approaches indeed lead to exactly the same 
AT. This underlines the validity of the indirect approach 
to obtain AT, which is often preferred in experimental 
studies, as long as the system does not actually cross the 
first order phase transition, i.e. for temperatures and elec¬ 
tric field strengths above the critical point We point 
out, though, that errors can be introduced due to imper¬ 
fect fits to the statistical measurements, in particular for 
small fields close to the phase transition, where polariza¬ 
tion and specific heat vary strongly. Another source of 
inaccuracies is the finite sampling of quantities as func¬ 
tion of temperature and field, which leads to numerical 
errors when integrating Eq. ©■ 

The isothermal entropy change calculated using direct 
(Eq. ©) and indirect (Eq. ([3])) methods is shown in 
Fig. Mb). We note that in principle the specific heat 
in Eq. (|2]) is the total specific heat, including electronic, 
ferroelectric, and all other lattice contributions, while in 
our treatment using the effective Hamiltonian, only the 
contributions from the ferroelectric soft mode variables 
are taken into account. However, since we use exactly the 
same degrees of freedom to also obtain the temperature- 
and field-dependent electric polarization for evaluating 
Eq. ([3]) , we obtain a consistent description within the ef¬ 
fective Hamiltonian, and both equations should in prin- 















ciple lead to the same value of AS. Nevertheless, in 
contrast to the adiabatic temperature changes, the cal¬ 
culated isothermal entropy changes do not need to be 
rescaled for the missing degrees of freedom in order to 
compare with experimental measurements. The entropy 
change in Eq. © depends only on differences in the spe¬ 
cific heat for different electric fields, and one can expect, 
at least to a good approximation, that only the soft mode 
variables will give significant contributions to this elec¬ 
tric field dependence. Consequently, also Eq. m does 
not contain any quantities depending explicitly on the 
missing degrees of freedom. The effect of the electrons 
and of other structural degrees of freedom on the tem¬ 
perature and field dependence of the polarization is, to 
a good approximation, implicitly taken into account by 
the soft mode variable. 

The entropy change shown in Fig. ©b) corresponds 
to a change of the electric field from 75kV/cm to 
300kV/cm, i.e. “switching on”, and therefore a nega¬ 
tive AS is obtained. Again, we exclude the region of 
small electric fields close to the phase transition for the 
comparison of direct and indirect methods. The lower 
bound for the temperature integration in Eq. 0 is cho¬ 
sen as Ti = 200 K. This temperature is above the sec¬ 
ond phase transition from the tetragonal ferroelectric to 
the orthorhombic ferroelectric phase in BaTiC >3 (with the 
chosen parameterization of the effective Hamiltonian). 
By definition, the “direct” AS calculated from Eq. 0 is 
zero for T = Ti, while the “indirect” A5 obtained from 
Eq. © has a finite value at T\ = 200 K. Since accord¬ 
ing to Fig. [1] the calculated specific heat at 200 K shows 
only negligible field dependence, the finite value of AS 
at this temperature is related to electric-field dependence 
of the specific heat at lower temperatures, most likely at 
the two ferroelectric-ferroelectric transitions (tetragonal- 
orthorhombic and orthorhombic-rhombohedral). For a 
better comparison between direct and indirect methods 
we therefore rigidly shift the AS curve obtained from 
the direct method such that it matches the AS value 
obtained from the indirect method at the lowest temper¬ 
ature Ti = 200 K (blue dotted line in Fig. ©. It can be 
seen that the shifted data agrees quite well with the data 
obtained from the indirect method. Small deviations can 
be observed close to the peak at around 350 K, which we 
ascribe to inaccuracies related to the smoothing/fitting 
of the specific heat and polarization data and to inte¬ 
gration errors due to finite temperature and electric field 
sampling. 

Interestingly, the obtained peak value of 
AS ~ 7 J-kg _1 K _1 is of the same order of magni¬ 
tude as the maximal value reported for BaTiC >3 single 
crystals measured in Ref. |29| (AS = 2.1 J-kg _1 K _1 ). 
However, the corresponding electric field intervals are 
completely different (75 to 300kV/cm in our calcu¬ 
lations, compared to 0 to 4kV/cm in Ref. [29h . and 
thus a meaningful quantitative comparison is not easily 
possible. 

Finally, in Fig. ©we compare the calculated adiabatic 


temperature change for different electric held intervals 
with the same width £; — £f| = 50kV/cm but differ¬ 
ent magnitude of S\ and £f. A first order phase tran¬ 
sition occurs only in the electric Held interval between 
£\ = 50kV/cm and £{ = OkV/cm. For this interval, 
we obtain a very narrow peak in AT at 310 K, with a 
maximum value of 15.8 K (corresponding to ~3.2K after 
scaling to the correct Nf). For larger applied Helds, the 
maximum AT value shifts to higher temperatures and 
the corresponding peak broadens. 

We note that from our specific heat calculations we 
can estimate a critical electric field of around 25 kV/cm. 
For larger fields, the polarization varies continuously with 
temperature, i.e. no first order phase transition occurs, 
and the temperature at which \dP/dT\ £ is maximal fol¬ 
lows the so-called “Widom line” (see e.g. Refs .Inland 1301. 
It can be seen that, while the largest AT is observed in 
the field interval containing the first order phase tran¬ 
sition (i.e. Si = 50kV/cm and £{ = OkV/cm), the field 
intervals corresponding to larger S\/S{ also give sizable 
contributions to the EC effect. We therefore conclude 
that while the vicinity to the first order transition is im¬ 
portant to obtain large changes of polarization with tem¬ 
perature and electric field, and thus large EC effect, the 
contribution of the transition itself is not essential to ob¬ 
tain large EC response. We note that similar conclusions 
have been reached in Ref. [30|, based on MD simulations 
for LiNbC> 3 . 

Further, a comparison between AT obtained using the 
direct and indirect methods for the different field inter¬ 
vals again shows a very good agreement for the field in¬ 
tervals corresponding to larger magnitude of £) and £f , 
where the variation with temperature and electric field 
is less strong. Clear discrepancies can be seen for the 
interval between 100 and 50 kV/cm, where AT is rather 
sharply peaked. These discrepancies result from imper¬ 
fect smoothing/fitting as well as from numerical integra¬ 
tion errors, as already discussed above. 


IV. SUMMARY AND CONCLUSIONS 

In summary, we have presented a computational study 
of the EC effect in BaTiC >3 using MD for a first principles- 
based effective Hamiltonian. We have compared the EC 
temperature change calculated using direct and indirect 
methods for bulk BaTiOs, thereby paying particular at¬ 
tention to the internal consistency of the method. In par¬ 
ticular, the temperature and electric field-dependent spe¬ 
cific heat has been calculated within the same framework 
as the temperature- and field-dependent electric polar¬ 
ization (required for the indirect determinations of AT), 
and the same framework has also been used for the direct 
calculation of AT using microcanonical MD. 

We have demonstrated that the direct and indirect de¬ 
termination of the adiabatic temperature change leads 
to identical results provided that the field and tempera¬ 
ture region very close to the first order transition, where 
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the indirect method is not applicable, is excluded. We 
note that the applicability of Maxwell’s relation, Eq. m, 
which underlies the indirect determination of the EC 
temperature change, has been critically discussed for sys¬ 
tems close to a first order phase transition (see e.g. the 
discussion in Refs. 0, [g, and 0). Our results clearly 
demonstrate the validity of this relation as long as the 
first order transition is not crossed. Directly at the first 
order transition, the specific heat diverges and dP/dT is 
not defined, and thus the indirect method is not appli¬ 
cable. Very close to the transition, errors can arise due 
to inaccurate fits, the use of a temperature- and field- 
independent specific heat, and due to finite temperature 
and field sampling of the integral in Eq. 0- 

Furthermore, we have demonstrated the importance of 
maintaining thermal equilibrium during the MD simu¬ 
lations for the direct calculation for AT, and we have 
shown that in the present case a ramping rate for the elec¬ 
tric field of 0.002 kV/cm/fs is sufficiently slow to ensure 
reversibility. We note, however, that, due to the neglect 
of the less important degrees of freedom in the effective 
Hamiltonian, this is not necessarily representative for the 
intrinsic relaxation time of BaTiOa. 

In addition, we have (to the best of our knowledge for 
the first time) used the effective Hamiltonian approach 
to calculate the isothermal EC entropy change. Similarly 
to the case of the adiabatic temperature change, we have 
found good agreement between (quasi-) direct, i.e. via the 
specific heat, and indirect determination of AS. While 
our calculated values are quantitatively of similar mag¬ 
nitude as available experimental data, further studies for 
different electric field strengths, possibly also consider¬ 
ing the contribution stemming from the latent heat of 
the first order phase transition, are necessary to obtain 
a more quantitative comparison between calculated and 
measured data. 

The observation that the largest EC temperature 


change occurs in the field interval containing the first 
order ferroelectric transition (£, = 50kV/cm and £f = 
OkV/cm, see Fig. [6j is in agreement to the giant 
caloric temperature changes found at other coupled 
ferroic-structural transitions, e.g. in magnetic Heusler al¬ 
loys In all these cases, small external fields are able 

to induce large adiabatic temperature changes, which, 
however, are restricted to only a narrow temperature in¬ 
terval. Unfortunately, in many cases the thermal hystere¬ 
sis of the transition leads to a significant reduction of the 
achievable reversible temperature changes under cycling 
of the fields, see e.g. Ref. [32]. In this respect, the EC ef¬ 
fect in BaTiC> 3 , with its rather low critical field strength, 
has an important advantage for cooling applications com¬ 
pared to the well-established magnetocaloric Heusler al¬ 
loys. Beyond the critical field strength, the thermal hys¬ 
teresis vanishes, leading to a fully reversible EC effect 
(see also the experimental results in Ref. (30). Further¬ 
more, the transition itself is not crucial for obtaining a 
large caloric response, as field intervals corresponding to 
larger £-J£i, i.e. above the critical field strength, also give 
sizable contributions to the EC effect, see Fig. [G] This is 
accompanied by a broadening of the AT peak with tem¬ 
perature, which is also advantageous for applications. In 
contrast, the strong first order character of the magneto- 
structural phase transition in magnetic Heusler alloys is 
conserved even for giant fields up to 40 T, without a re¬ 
duction of the thermal hysteresis^ 
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